#ifndef _GEODESICS_H_
#define _GEODESICS_H_
//This file implement the modified Euler Method to compute the geodesics after applying FMM on the manifold
//#include <boost/numeric/mtl/config.hpp>
#include <boost/numeric/mtl/mtl.hpp>


typedef mtl::dense2D<double>  MATRIX;
typedef mtl::Collection<MATRIX>::value_type SCALAR;
typedef mtl::dense_vector<SCALAR>            VECTOR;

bool LU_solver (const MATRIX & A, const VECTOR & b, VECTOR & solution);
bool three_points_interpolation(Mesh::Point v[3], Mesh::Point & intersection, int next_face);
bool six_points_interpolation (Mesh::Point & start, Mesh::Point v[], double value[], double step, Mesh::Point & intsec, int & next_face);
void draw_geodesics (Mesh & mesh, const Mesh::Point & start, Mesh::FaceHandle start_face,  const Mesh::Point & end, 
					 const std::vector<Mesh::FaceHandle> & end_faces, double step, std::vector<Mesh::Point> & path);
void generate_merged_mesh_geodesics(Mesh & mesh, std::vector<Mesh::Point> & path);
#endif
